*! by Hung-Jen Wang
*!    Department of Economics, National Taiwan University
*!    wangh@ntu.edu.tw


** This Stata program estimates parameters of a SF frontier model and performs
** the sine and cosine tests on the distribution assumptions of the composed error.
** The null hypothesis can be either N-HN or N-Exponential distribution.
** It also optionally computes the JLMS and BC (in)efficiency index.
** The estimation is based on the second and third moment restrictions and the test is based on
** the characteristic function restrictions.
** Please see 00Readme.txt for more information.
**
** Usage examples:
**        sfmtest y, prod frontier(x1 x2 x3) model(h) omega(1)
**        sfmtest y, cost frontier(x1 x2 x3) model(e) omega(1)
**        sfmtest y, prod frontier(x1 x2 x3) model(h) omega(1.5) bc(myindex1) jlms(myindex2)
**
** Citation:
**   Chen, Y.-T., and Wang, H.-J. (2012) "Centered-Residuals-Based Moment Estimator and
**   Test for Stochastic Frontier Models", Econometric Reviews 31(6), pp.625-653.
**
** program contact: Hung-Jen Wang, wangh@ntu.edu.tw
**                  http://homepage.ntu.edu.tw/~wangh
**                  2012/04/14


capture program drop sfmtest
program define sfmtest

version 10.1

syntax varlist [if] [in], FRONTIER(string) UDISTribution(string) [ PRODuction COST OMEGA(string) JLMS(string) BC(string) UCOE(string) VCOE(string) USTD(string) VSTD(string) Cosine(string) Sine(string) NOREPOrt ]
marksample touse, nov

preserve /* keep the current copy of the dataset, so that in case of -cost-, do not need to switch sign back when estimation is done or interrupted */

tempvar Resid CResid dm2 dm3 Grad2 Grad3 TFc TFs Gradc Grads xic xis xic2 xis2 xics2
tempvar Grad22 Grad33 Grad23
tempvar useobs

tempname m2 m3 lambda2 lambda3 s2u s2v w Ecoswv EDcoswv mu Dmu Ecoswu Edcoswu
tempname t erfiw Esinwu EDsinwu hc Mc Dhc hs Ms sT
tempname Dhs h11 h12 h21 h22 dc1 dc2 ds1 ds2 mTFc mTFs xic2sum Testc xis2sum Tests inv11 inv12 inv22 Testcs
tempname mGrad22 mGrad33 mGrad23

/* --- determine whether it's a cost or a production frontier --- */

if  (("`cost'"~="") & ("`production'"~="")) {
     di in red "Need to specify one and only one of the -cost- and -production- option."
     error 200
}

***********************
  /* --- flip data if estimate cost frontier ------- */

if "`cost'" ~= "" { /* estimate a cost frontier */

  unab fronvar: `frontier'

   foreach X of local fronvar  {
     quie replace `X' = -`X'
   }
     quie replace `varlist' = -`varlist'
}

****************************

quie regress `varlist' `frontier' if `touse'
quie gen `useobs' = cond(e(sample), 1, 0) /* it is based on touse and the regression */
quie predict double `Resid' if `useobs', resid

quie sum `varlist' if `useobs'
scalar `sT' = r(N)


quie summarize `Resid' if `useobs'
quie gen double `CResid' = `Resid'- r(mean)  if `useobs' /* center residual */

    gen double `dm2' = (`CResid')^2  if `useobs'
    quie summarize `dm2'  if `useobs'
    scalar `m2' = r(mean)
    gen double `dm3' = (`CResid')^3  if `useobs'
    quie summarize `dm3'  if `useobs'
    scalar `m3' = r(mean)

  if "`udistribution'" == "h" { /* null is half-normal */
    scalar `lambda2' = 1-2/(_pi)
    scalar `lambda3' = (1-4/(_pi))*sqrt((2/(_pi)))
  }
  if "`udistribution'" == "e" { /* null is exponential */
    scalar `lambda2' = 1
    scalar `lambda3' = -2
  }

    scalar `s2u' = ((`m3'/`lambda3')^2)^(1/3) /* estimated sigma_u^2 */
    scalar `s2v' = `m2'-`lambda2'*`s2u'   /* estimated sigma_v^2 */


*** prepare for later report ****

quie regress `varlist' `frontier'  if `useobs'


 if "`udistribution'" == "h" {
   quie nlcom _b[_cons] + sqrt(`s2u')*sqrt(2/_pi)
 }
 if "`udistribution'" == "e" {
   quie nlcom _b[_cons] + sqrt(`s2u')
 }

if "`cost'" == "" {
  local sign = -1
}
else {
  local sign = 1
}

 mat b0_olsM = r(b)
 mat b1_olsM = r(V)

 local b0_ols = `sign'*(b0_olsM[1,1])
 local b1_ols = sqrt(b1_olsM[1,1])


*************

local myucoe = `s2u'
local myvcoe = `s2v'

if "`ucoe'" ~= "" {
   scalar `ucoe' = `s2u'
 }
if "`vcoe'" ~= "" {
   scalar `vcoe' = `s2v'
 }

    gen double `Grad2' =(`CResid')^2-(`s2v'+`lambda2'*`s2u')  if `useobs'
    gen double `Grad3' =(`CResid')^3-`lambda3'*((`s2u')^1.5)-3*(`s2v'+`lambda2'*`s2u')*(`CResid')  if `useobs'
    mat HessE = (-1, -`lambda2' \ 0 , -1.5*`lambda3'*sqrt(`s2u'))


/* --- Var-Cov matrix for the estiamtes ---- */


quie gen double `Grad22' = (`Grad2')^2  if `useobs'
quie gen double `Grad33' = (`Grad3')^2  if `useobs'
quie gen double `Grad23' = (`Grad2'*`Grad3')  if `useobs'

quie sum `Grad22'  if `useobs'
 scalar `mGrad22' = r(mean)
quie sum `Grad33'  if `useobs'
 scalar `mGrad33' = r(mean)
quie sum `Grad23'  if `useobs'
 scalar `mGrad23' = r(mean)

mat mMat = (`mGrad22', `mGrad23' \ `mGrad23', `mGrad33')
mat vcovM  = (1/(`sT'))*inv(HessE)*mMat*inv(HessE)'

local myvstd = sqrt(vcovM[1,1])
local myustd = sqrt(vcovM[2,2])

if "`vstd'" ~= "" {
  scalar `vstd' = sqrt(vcovM[1,1])
}
if "`ustd'" ~= "" {
  scalar `ustd' = sqrt(vcovM[2,2])
}



if "`omega'" == "" { /* just so the program runs */
 scalar `w' = 1
 local omega = 1
 }
else {
 scalar `w' = `omega'
}


scalar `Ecoswv' = exp(-0.5*((`w')^2)*`s2v')
scalar `EDcoswv'= -0.5*((`w')^2)*`Ecoswv'


if "`udistribution'" == "h" { /* The normal-half normal model */
  scalar `mu' = (sqrt(`s2u'))*sqrt(2/(_pi))
  scalar `Dmu'= (sqrt(2/(_pi)))/(2*sqrt(`s2u'))
  scalar `Ecoswu' = exp(-0.5*((`w')^2)*`s2u')
  scalar `Edcoswu' = -0.5*((`w')^2)*`Ecoswu'
  scalar `t' = `w'*sqrt(`s2u'/2)
  myerfi, upper(scalar(`t')) generate(`erfiw')
  scalar `Esinwu' = `Ecoswu'*`erfiw'
  scalar `EDsinwu' = `Edcoswu'*`erfiw'+`w'/sqrt(2*(_pi)*`s2u')
}


if "`udistribution'" == "e" {/* The normal-exponential model */
  scalar `mu' = sqrt(`s2u')
  scalar `Dmu'= 1/(2*(sqrt(`s2u')))
  scalar `Ecoswu' = 1/(1+`s2u'*((`w')^2))
  scalar `Edcoswu' = -((`w')^2)/((1+`s2u'*((`w')^2))^2)
  scalar `Esinwu' = (`w'*(sqrt(`s2u')))/(1+`s2u'*((`w')^2))
  scalar `EDsinwu' = `w'/(2*(sqrt(`s2u'))*(1+`s2u'*((`w')^2)))-(((`w')^3)*(sqrt(`s2u')))/((1+`s2u'*((`w')^2))^2)
}




scalar `hc' = `Ecoswu'*cos(`w'*`mu')+`Esinwu'*sin(`w'*`mu')
scalar `Mc'=`Ecoswv'*`hc'
scalar `Dhc'=`Edcoswu'*cos(`w'*`mu')+`EDsinwu'*sin(`w'*`mu')-`Ecoswu'*sin(`w'*`mu')*`w'*`Dmu'+`Esinwu'*cos(`w'*`mu')*`w'*`Dmu'
mat Dc = (-`EDcoswv'*`hc' \ -`Ecoswv'*`Dhc')

scalar `hs' = `Ecoswu'*sin(`w'*`mu')-`Esinwu'*cos(`w'*`mu')
scalar `Ms' = `Ecoswv'*`hs'
scalar `Dhs' = `Edcoswu'*sin(`w'*`mu')-`EDsinwu'*cos(`w'*`mu')+`Ecoswu'*cos(`w'*`mu')*`w'*`Dmu'+`Esinwu'*sin(`w'*`mu')*`w'*`Dmu'
mat Ds = (-`EDcoswv'*`hs' \ -`Ecoswv'*`Dhs' )
gen double `TFc'=cos(`w'*(`CResid'))-`Mc'  if `useobs'
gen double `TFs'=sin(`w'*(`CResid'))-`Ms' if `useobs'
gen double `Gradc'=`TFc'+`Ms'*`w'*(`CResid') if `useobs'
gen double `Grads'=`TFs'-`Mc'*`w'*(`CResid') if `useobs'


/* Tests for the model. */


mat invHess = inv(HessE')
scalar `h11' = invHess[1,1]
scalar `h12' = invHess[1,2]
scalar `h21' = invHess[2,1]
scalar `h22' = invHess[2,2]


scalar `dc1' = Dc[1,1]
scalar `dc2' = Dc[2,1]
scalar `ds1' = Ds[1,1]
scalar `ds2' = Ds[2,1]

gen double `xic' = `Gradc' - (`dc1'*(`Grad2'*`h11'+`Grad3'*`h21')+`dc2'*(`Grad2'*`h12' + `Grad3'*`h22'))  if `useobs'
gen double `xis' = `Grads' - (`ds1'*(`Grad2'*`h11'+`Grad3'*`h21')+`ds2'*(`Grad2'*`h12' + `Grad3'*`h22')) if `useobs'

quie summarize `TFc' if `useobs'
scalar `mTFc' = r(mean)


quie summarize `TFs' if `useobs'
scalar `mTFs' = r(mean)

mat mTFcs = (`mTFc' \ `mTFs')


* mat CTestcNH = sT*mTFcNH*inv(xicNH'*xicNH/sT)*mTFcNH
quie gen double `xic2' = (`xic')^2 if `useobs'
quie sum `xic2' if `useobs'
scalar `xic2sum' =1/r(mean)
local mycosine = `sT'*`mTFc'*`xic2sum'*`mTFc'

if "`cosine'" ~= "" {
  scalar `cosine' = `mycosine'
}

* mat CTestsNH = sT*mTFsNH*inv(xisNH'*xisNH/sT)*mTFsNH
quie gen double `xis2' = (`xis')^2 if `useobs'
quie sum `xis2' if `useobs'
scalar `xis2sum' =1/r(mean)
local mysine = `sT'*`mTFs'*`xis2sum'*`mTFs'

if "`sine'" ~= "" {
 scalar `sine' = `mysine'
}

* mat CTestcsNH = sT*mTFcsNH'*inv(xicsNH'*xicsNH/sT)*mTFcsNH
quie gen double `xics2' = `xic'*`xis' if `useobs'
quie sum `xic2' if `useobs'
scalar `inv11' = r(sum)
quie sum `xis2' if `useobs'
scalar `inv22' = r(sum)
quie sum `xics2' if `useobs'
scalar `inv12' = r(sum)
mat invmat = (`inv11', `inv12' \ `inv12', `inv22')/`sT'

mat CTestcs = `sT'*mTFcs'*inv(invmat)*mTFcs

scalar `Testcs' = CTestcs[1,1]


***** JLMS *************************************


if ("`jlms'" ~= "") | ("`bc'" ~= "") {

restore  /* restore the data!! */

tempvar  myx mustar sig2star rat1 myeps useobs

quie regress `varlist' `frontier' if `touse'
quie gen `useobs' = cond(e(sample), 1, 0) /* it is based on touse and the regression */
quie predict double `myx' if `useobs', xb


 if "`udistribution'" == "h" {
   quie replace `myx' = `myx' + sqrt(`s2u')*sqrt(2/_pi) if `useobs'
   quie gen double `myeps' = `varlist' - `myx' if `useobs'
   if "`cost'" ~= "" { /* meaning a cost frontier */
     quie replace `myeps' = -1*`myeps' if `useobs'
   }
   quie gen double `mustar'  = -(`myucoe')*(`myeps')/(`myucoe' + `myvcoe') if `useobs'
   quie gen double `sig2star' = (`myucoe')*(`myvcoe')/(`myucoe' + `myvcoe') if `useobs'
   quie gen double `rat1' = (`mustar')/(sqrt(`sig2star')) if `useobs'
   if "`jlms'" ~= "" {
     quie gen double `jlms' = sqrt(`sig2star')*normalden(`rat1')/normal(`rat1') + (`mustar') if `useobs'
     label var `jlms' "conditional E(u|e)"
   }
   if "`bc'" ~= "" {
     quie gen double `bc' = exp(-`mustar' + 0.5*`sig2star')*normal(`rat1' - sqrt(`sig2star'))/normal(`rat1')
     label var `bc' "conditional E(exp(-u)|e)"
   }
 }
 if "`udistribution'" == "e" {
   quie replace `myx' = `myx' + sqrt(`s2u') if `useobs'
   quie gen double `myeps' = `varlist' - `myx' if `useobs'
   if "`cost'" ~= "" { /* meaning a cost frontier */
     quie replace `myeps' = -1*`myeps' if `useobs'
   }
   quie gen double `mustar' = -`myeps' - (`myvcoe')/(sqrt(`myucoe')) if `useobs'
   quie gen double `rat1' = (`mustar')/sqrt(`myvcoe') if `useobs'
   if "`jlms'" ~= "" {
     quie gen double `jlms' = sqrt(`myvcoe')*normalden(`rat1')/normal(`rat1') + `mustar' if `useobs'
     label var `jlms' "conditional E(u|e)"
   }
   if "`bc'" ~= "" {
     quie gen double `bc' = exp(-`mustar' + 0.5*(`myvcoe'))*normal(`rat1' - sqrt(`myvcoe'))/normal(`rat1') if `useobs'
     label var `bc' "conditional E(exp(-u)|e)"
   }
 }

}


***** critical values **************************

local crit1 = invchi2(1,1-0.01)
local crit2 = invchi2(1,1-0.05)
local crit3 = invchi2(1,1-0.10)


***********************************************

if "`noreport'"=="" { /* report results */

di " "
 di " "
if "`udistribution'" == "h" {
 di in yel "***********************************************************************"
 di in yel "**** Moment Estimates and Tests of the Normal - Half Normal Model ****"
 di in yel "***********************************************************************"
 }
else if "`udistribution'" == "e" {
 di in yel "***********************************************************************"
 di in yel "**** Moment Estimates and Tests of the Normal - Exponential Model ****"
 di in yel "***********************************************************************"
 }

 di in gre " "
 di "     intercept = " in yel `b0_ols' in gre ", std.err. = " in yel `b1_ols' in gre "."
 di "sigma-v-square = " in yel `myvcoe' in gre ", std.err. = " in yel `myvstd' in gre "."
 di "sigma-u-square = " in yel `myucoe' in gre ", std.err. = " in yel `myustd' in gre "."
 di " "
 di "Coefficients of other frontier variables can be obtained from OLS."



di " "
di " "
di in yel "  **** Sine and Cosine Tests of the Model's Composed Error (v-u or v+u) ***"
di in gre " "
if "`udistribution'" == "h" {
  di "Null: v is normal AND u is " in yel " half-normal" in gre "."
}
else if "`udistribution'" == "e" {
  di "Null: v is normal AND u is " in yel "exponential" in gre "."
}
di " "
di "With omega (w) = " in yel `omega' in gre ", the sine test statistic = " in yel `mysine' in gre ", and"
di in gre "                     the cosine test statistic = " in yel `mycosine' in gre "."
di " "
di "The 1%, 5%, and 10% critical values of the chi-square distribution with"
di "1 d.o.f. are, respectively, " in yel `crit1' in gre ", " in yel `crit2' in gre ", and " in yel `crit3' in gre "."
di " "

if "`jlms'"~="" | "`bc'" ~= "" {
  di " "
  di " "
  di in yel "  **** The JLMS and/or BC (in)efficiency index ***"
  di in gre " "
  if "`jlms'" ~= "" {
    di "The JLMS index ( E[u|epsilon] ) is saved to " in yel "`jlms'" in gre "."
  }
  if "`bc'" ~= "" {
    di "The BC index ( E[exp(-u)|epsilon] ) is saved to " in yel "`bc'" in gre "."
  }
}

}


end
